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The authors present a study of the non equihbrium statistical properties of a one dimensional 
hard-rod fluid dissipating energy via inelastic collisions and subject to the action of a Gaussian heat 
bath, simulating an external driving mechanism. They show that the description of the fluid based 
on the one-particle phase-space reduced distribution function, in principle necessary because of the 
presence of velocity dependent coUisional dissipation, can be contracted to a simpler description 
in configurational space. Indeed, by means of a multiple-time scale method the authors derive a 
self-consistent governing equation for the particle density distribution function. This equation is 
similar to the dynamic density functional equation employed in the study of colloids, but contains 
additional terms taking into account the inelastic nature of the fluid. Such terms cannot be derived 
from a Liapunov generating functional and contribute not only to the relaxational properties, but 
also to the non equilibrium steady state properties. A validation of the theory against molecular 
dynamics simulations is presented in a series of cases, and good agreement is found. 

PACS numbers: 02.50.Ey, 05.20.Dd, 81.05.Rm 



I. INTRODUCTION 



Granular fluids (GFs) represent one of the current paradigms of open non-equilibrium systems and, for this reason, 
in the last two decades have been the subject of a huge amount of experimental, numerical and theoretical stud- 
jggij2j3^i5j6j7^ GFs can be conveniently modeled as assemblies of macroscopic particles, experiencing instantaneous 
binary collisions during which a fraction of the kinetic energy is dissipated, i.e., transferred into internal degrees of 
freedom. Under the action of a vigorous external driving force, GFs may appear similar to ordinary molecular fluids, 
but crucial differences remain because inelasticity leads to the appearance of a series of peculiar behaviors, such as 
clustering, non-Gaussian velocity distribution, and velocity correlations. These phenomena have no counterparts in 
molecular fluids and render the study of GFs difficult but particularly fascinating. 

In spatially uniform systems, relations have been obtained between static average quantities such as density, kinetic 
temperature and pressure, which may be regarded as the analogue of the equation of state. In addition, a granular 
hydrodynamics has been developed which, due to the inelasticity of collisions, differs nontrivially from standard 
hydrodynamics. The majority of these studies focus on large scale properties of the fluid. However, in strongly 
inhomogeneous systems, the connection between the microscale typical of the particles and the macroscale is still 
incomplete. Recently, some authors^ have proposed phenomenological theories, based on local mass and momentum 
conservation laws, incorporating nonideal gas effects via an effective free energy functional suitably designed to describe 
the spontaneous formation of loosely and densely packed regions. This approach sounds very appealing because the 
free energy density functional^, besides being a method computationally simple and physically clear, has proven to 
be a useful tool in the theory of nonuniform fluids with applications to interfacial and freezing phenomena. The 
basic assumption of all density functional theories is that the thermodynamic potential of a nonuniform system may 
be approximated knowing the structural and thermodynamic properties of the corresponding uniform system. Two 
questions are in order before proceeding to generalize the density functional theory (DFT) to granular materials: Does 
the same method offer any new insight in this new area? How far meaningful concepts for standard molecular fluids, 
such as free energy and chemical potential, can be extended to systems which are not at thermodynamic equilibrium? 
The answer to the second question seems to be desolately negative and therefore in order to construct a theory of 
nonuniform GF, alternative techniques not involving free energy functional derivations have to be developed. A step 
towards this new direction has been recently made and an equation of evolution for the particle density was obtained 
in the case of nonequilibrium colloidal systems^''. Interestingly, such an equation is very similar to the equation 
considered within the dynamic density functional method^!. It differs from the latter only for the presence of some 
higher order corrections. The derivation of Ref.^° employs standard tools of kinetic theory and in particular the 
revised Enskog equatio n^^i^^ . In this paper we generalize such an approach to fluids of inelastic hard core particles 
subject to a stochastic forcing. 
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We consider a one dimensional model of granular fluid which is simple enough as to lend itself to analytic work, but is 
endowed with a sufficient complexity as to display inhomogeneous behavio r-'^'*i-'^^'^^i^'''i'^^i-'^^'^°'^-'^i^^'^'^i^'^ . One dimensional 
models may play a useful role since they can be employed to test approximations of more general applicability and 
allow us to link easily the structural properties to the dynamical behavior. Moreover, at equilibrium the qualitative 
differences between one dimension and higher dimensions appear in the development of long-range ordering or any 
kind of phase transitions, but not in the short-range packing structure. In the collision kinetics, the qualitative 
difference between one dimension and higher dimensions appears in purely inertial and elastic systems through the 
conservation of the velocity distributions despite the collisions. The role of the bath, and also the inelastic collisions, 
kills the peculiarity of the one dimensional case. A basic feature of this work is the assumption of a uniform thermostat 
to describe the external energy supply. The balance between thermostatting and dissipation mechanism gives rise 
to noncquilibrium steady states which are achieved without fine tuning of the model parameters. Few kinds of 
thermostats have been employed in the literature, namely, the white noise thermostat - the Gaussian thermostats^ 
and the Langevin thermosta t^^i^^'^^ which includes both the white noise term and the friction force proportional to 
the velocity of the particles^. The present study is based on the Langevin thermostat because it can be easily realized 
in numerical experiments and lends itself to a great deal of analytical work. 

The paper is organized as follows. In Sec. [TTl we introduce the equations describing the dynamics of the stochastically 
driven inelastic hard-rod fluid model. We start from the stochastic equations for the trajectories of each particle and 
discuss how these can be reduced under some suitable approximations to the Fokker-Planck-Boltzmann equation for 
the single-particle phase-space distribution. Next, in order to make analytical progress, we separate the velocity and 
the spatial dependence of the distribution function and obtain an infinite hierarchy of coupled integro-differential 
equations. In Sec. IIIII we analyze the steady state uniform properties of the system. In Sec. IIVI we introduce the 
evolution equation for the density of the system, which is obtained in Appendix A by applying the multiple time 
scale method, and in Sec.|V]we perform a series of numerical tests of our theory using a few selected examples. The 
tests compare the two levels of description: the results obtained studying the single particle trajectories with those 
obtained from the density equation. Finally, conclusions are drawn in Sec. IVII 



II. MODEL 



Let us consider a one dimensional fluid consisting of N identical inelastic hard rods of mass m, length a, coefficient 
of restitution a, positions Xi, and velocities Vi, with i = 1,...,N. For the sake of generality we also consider an 
arbitrary external force, fext{x). When the separation of particles i and ? -I- 1 is cr a binary inelastic collision occurs. 
The collision conserves the total momentum and is described by the linear transformation: 

v'l = vi + ^-^(^'2 ~vi) (1) 

' l + a 
^2 = "^^2 {V2 - vi) 

connecting the precoUisional velocities (primed symbols) and the postcoUisional velocities (unprimed symbols). Since 
in a single collision the amount of total kinetic energy, 

AE^-^{l-a')iv[-v',)\ (2) 

is dissipated, without energy injection the particles would come to rest. On the contrary, a steady state regime can 
be reached if the energy loss through collisions is balanced by an energy injection that we assume to be realized 
by the combination of a friction force —mjVi and a stochastic force £,i{t). These two forces represent the Langevin 
thermostat. The complete dynamics can be represented by the following 2N coupled stochastic differential equations: 

^ (3) 

= _rn-fv, + f,^t{x^) + fr" + m (4) 

where indicates symbolically the resultant of the impulsive forces acting on particle i in possible hard core 

collisions against other particles. The stochastic force ^i(i) has zero average {(,i{t)) and white noise correlation 

Mt)CAt')) = 2^mToS,,S{t - t') , (5) 

The amplitude Tq is the "heat-bath temperature" and (•) indicates the average over a statistical ensemble of realizations 
of the noise. A statistical description of the system in terms of the one-particle phase space distribution /^^^ (x, d, t). 
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giving the number of particles in the volume element {x, x + dx, f, f + dw), can be worked out by taking the average 
over all realizations of the stochastic noise (see. RefJ^). Moreover, the distribution /'^^^(x,v,t) evolves according to 
the governing equation 



^ fext(x) d 

dx 



dv 



d_ 

dv 



m dv^ 



f^^\x,v,t) + k{x,v,t). 



(6) 



In the left hand side, the term between the square brackets describes the free streaming of the particles subject 
to the external force fe{x), the first term in the right hand side is the one-particle Fokker-Planck collision term 
representing the interaction with the heat bath, while k(x, v, t) describes the collisions among the particles. We treat 
these interactions within the revised Enskog theory (RET), developed by Ernst and van Beijerer>i2,. The RET for 
elastic collisions is accurate over the entire fluid range and describes the crystal phase too. It has been generalized to 
the inelastic regime and used to derive transport coefficients^!^. We write the RET collision operator as: 



k{x,v,t)= ^ dv2Q{svi2){svi2) 

LI J 



(7) 



s=±l 



r 1 



g2{x,x- sa\p)f'^^\x,v[,t)f^'^'>{x - sa, v'2,t)- g2{x,x + sa\p)f^^\x,vi,t)f^^'>{x + sa, V2,t) 



Notice that at variance with the elastic case, a quadratic factor in the gain term, specific to granular gases, appears. 
One power is the consequence of the Jacobian dv[dv2 — ^dvidv2 and the second power stems from the reflection law 
vi2 = —av'i2- The sum over s = ±1 is the analogue in d = 1 of the integration over the d-dimensional surface of the 
hyper-sphere with radius a. The RET embodies spatial correlations through the hard-rod pair correlation function, 
g2ix, x±a;n) evaluated at contact. As a simplifying approximation, to obtain a theory at the level of the one-particle 
distribution, we take g2{x, x ± a\p) to be given by its equilibrium valued evaluated when the local density is p(x, t): 



g2ix±a\p) 



1 



1-Vix±i)' 



(8) 



The time and density dependence occurs entirely via the local packing fraction ri[x, t) = /^J^J^j^ dx' p{x' , t). Therefore, 
the collision operator Eq. ^ is approximated by an explicit non-local functional of the one-particle density distribu- 
tion, with the terms Z'^-* (x, f i, i)/*-^-* (x ± cr, U2, i), set by the collision distance, and the nonlocal density dependence 
through ri{x ± cr/2), to include the particle correlations. 

In the following we shall employ the non dime nsional set of variables which are obtained by measuring the velocities 
in units of the thermal velocity = y^kBTo/m and lengths in unit of cr, i.e. V = v/vt and X = x/cr. The remaining 
variables can be non-dimensionalized according to the transformations r = tvx/sigma), T = ^a/vx- F(X) = 
<^ fextix) /mv'^. Finally, the distribution function and the collision term are rescaled according to the transformations: 
P(X,y,T) = avTf^^Kx,v,t) and K{X,V,t) = a^k{x,v,t). 

Equation ^ can be cast in the following non dimensional form: 



idP{x,v,T) 
r" 



dr 



LfpP{X, V, t) - it/ Ap(x, r) - ^F{X, r)-^P{X, V. r) + ^K{X, V, r) (9) 



where we have introduced Fokker-Planck operator, Lpp by the equation: 



The eigenfunctions oi Lpp read explicitly 



P{X,V,t) 



i^.(^)-^(-l)^^exp(--^^) 



(10) 



(11) 



and correspond to discrete eigenvalues v ~ 0, —1, —2, .. . We separate the velocity from the spatial dependence by 
expanding, over the basis set Hi,{V), both the phase-space distribution 



P(X,t/,T)=^<i>,(X,T)i/,(l^) 



(12) 
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and the collision term 



T/, r) = ^ C.(X, T)iJ,(y). 



(13) 



As shown in Appendix A, the coefficients C,^{X,t) can be expressed as nonlocal products of the moments r) 
with coefBcients which are nonlocal functionals of the density distribution. Substituting Egs. p^ and (fT51) in Eq. ^ 
and using the orthogonality of the basis set H,y{V), we obtain a system of coupled equations: for the moments ^p{X, r) 
which can be written in compact form as: 



dr 



dX 



L(X,r) =0 



(14) 



with $_! = 0. 

We identify the moment $o(^,t") = J dVP{X,V,T) with the number density, $i(X,t) = / dVVP{X,V,T) with 
the momentum density and ^2{^,t) + ^o{X,t)/2 = 1/2 J dVV'^P{X,V,T) with the kinetic energy density. For 
v ^ 0,1,2, Eq. ([H]) encodes the balance equations for these moments, i.e. the hydrodynamic equations characterizing 
a viscous oncdimensional fluid^'*. Finally, we introduce a local kinetic temperature, often called granular temperature, 
via the definition: 



T{X,T) = {{V^)^{Vf) = l + 2 



^2(^,r) 



^i{X, t) 



(15) 



III. HOMOGENEOUS STEADY STATE PROPERTIES 



Before embarking upon the task of solving the evolution equations, we illustrate the peculiarity of the inelastic 
system by choosing the simplest case, namely, a time-independent spatially uniform system with F{X) = 0. We 
consider the global velocity distribution function '5(1^), and show that in the steady state it does not relax to the 
Maxwellian, as it would occur in the case of a molecular fluid. To this purpose let $i = $i, where $i are some 
constants. The value of the amplitude $2, through Eq. (|14p . can be expressed in terms of the miiform density $0: 



2r$2 = C2 = -^^-^g2[H + 3$o^2 + 7^2] : 

VTT 4 



(16) 



where the second equality follows from the definition of C2 and the expressions of appendix A. To first order in the 
inelasticity parameter (1 — a^) we obtain: $2 = — e$o/2, where e = [(1 — a^) / {y/7TT)]g2^o- The procedure can be 
carried on for values of v larger than 2 with the following result: 



<I>4 



4r 96 



(17) 



' 6r ^ 5760 



- 1 ^ (15 - 20^2 + 50a4 - 16a6) ~ 
- 8f ^« ^ 215040 



Hence, ^'(T^) can be written as 



$0 + [V^ - 1)$2 + {V^ - + 3)$4 



{V^ - 15y^ + 45^2 _ I5)$g + (ys _ 28^^ + 210^^ - 420^^ + 105)$5 



(18) 
(19) 

(20) 



Two remarks are in order: if we retain only the two leading terms in the expansion, the distribution function can be 
approximately rewritten as a Maxwellian, 



^{V) = ^o l-(l^'-l)- 



-^V^^ exp(-^) ^ 
/2^ ^27r(l - e) 



(21) 
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and we interpret e as the reduced temperature shift induced by the inelastic dissipation. Secondly, the expansion 
Eq. (pij) can be compared with an exact solution of Eq. valid when F{X) = and in the limit (1 — a) — > 0"^^, 
obtained by Benedetto et al.. These authors showed that a spatially uniform solution, 'i>p{V), of Eq. ([9]) is given 
implicitly by the following nonlinear integral equation: 



exp<^ 



(1 - a)g2 



6r 



duu^-^piu + V) - I duu^^p{u + V) I 



(22) 



where Z is the constant which ensures the correct normalization of the probability distribution function (PDF). 
Interestingly, such a distribution has high- velocity tails which decay as exp(— c|yp), whereas the central region of the 
distribution is approximately a Maxwellian. Clearly, the high-velocity tails cannot be well reproduced by the present 
expansion, which is applicable when F >> 1, but the kurtosis associated with Eq. (|2ip compares reasonably with the 
kurtosis computed from the distribution Pp{V), as shown in Fig. [T] 
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FIG. 1: Kurtosis of the velocity distribution as a function of a for po = 0.8 and F = 5. The dashed line represents the data of 
the present theory, whilst the continuous line refers to the results from formula (|22p by Pulvirenti and coworkers. In the inset 
we display the corresponding kinetic temperatures. 



IV. EVOLUTION EQUATION 

We shall consider, in the following, the nearly overdamped regime F > i36, 37,38, 39, 4o,4i,42 _ Since only the particle 
number is conserved, one expects that after a transient of duration of the order of 7~^, the momentum and the energy 
current become slaved by the density field. This remark allows us to simplify the task posed by the open hierarchy of 
Eqs. (HI]). In Ref.^°, we showed that, for a system undergoing perfectly elastic collisions, the problem can be treated 
conveniently by employing a multiple-time scale technique. As a result we found a reaction-diffusion self-consistent 
equation involving only the amplitude ^{){X^t). The evolution of all remaining partial amplitudes <I>^(X, t), [y > 1) 
could be deduced from the knowledge of ^o{X,t). Physically, the reason for such a complexity reduction can be 
attributed to the fact that the marginal velocity probability distribution attains its local equilibrium rapidly, in a 
time span of the order of 7~^, during which the one-particle density changes slowly. Indeed, the positional degrees of 
freedom reach an equilibrium distribution on a much slower time scale than the velocities. 

Since the method of solution follows closely the derivation of Refii^, we report the details of the present case in 
Appendix B and proceed to illustrate the resulting equation of evolution. We only recall that the method is based 
on a systematic expansion in powers of F~^ which takes into account the fact that in Eq. ^ the time derivative is 
multiplied by the small parameter F^^. For such a reason a multiple time-scale method has to be applied. We also 
need to introduce the following expansions of the moments and of the collision integrals in inverse powers of F: 

$,(X,r)=^i^(/.„.(X,T) (23) 

n 

and 

a(X,T)=^^c„,(X,r) (24) 
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The key result of the analysis contained in Appendix B is the following equation for the density amplitude 



del), 



'00 



dr 



^dx{[dx - FiX)]MX,T) - coiiX,T) - ^cn(X,T) + ^dxCo2{X,T)} . 



(25) 



Equation (|25p is the fundamental equation of this work and constitutes a closed expression, once the collisional terms 
Csi, and the amplitudes (jisf are specified in terms of the scaled density ^oo ''') • -A-S we will show below, the density 
field (j)oQ{X,T) fully characterizes the state of the system and slaves the remaining hydrodynamic fields. Indeed, the 
amplitudes of the Hi{V) and H2{V) components are completely determined from the knowledge of (j)oo{X,T). Such 
a complexity reduction occurs because the density is the only conserved field in our thermostatted model. 
Using Eq. (jB4p we obtain at order r~-^, the following: 



coi{X,t) 



(1 + a) 



(f>oo{X, t) 52 (X, X + l)0oo(^ + 1, r) - g2(^, ^ - l)0oo(^ - 1, r) 



and the following at order V ^: 



C02(^,t) = -- 



000 {X, t) 52 (X, X + l)(/.oo (X + 1 , r) + 52 (X, X - l)0oo (X - 1 , r) 



(26) 



(27) 



and 



Cii(X, r) 



(1 + a) 



000 (X,t) 52(X,X + 1) 



(1 



-0ii(X,t) 52(X,X + 1) 



hi{X- 
^00 (X - 



1,t) 



1) 



1,t) 



1) 



2 

(1 + a) 



.92(X,X 

.92(X,X 

ha{X, t) [g2{X, X + l)0i2(X + 1, r) - 52(X, X - 1 
h2{X, t) L(X, X + l)0oo(X + 1, r) - 52(X, X - 1 



ii{X- 
m{X 



1,t) 



1,t) 



/.i2(X-l,r) 
/.oo(X- l,r) 



where 



and 



011 (X, r) = - [ax - F(X)]0oo(X, r) + coi(X, r) 



1 



012(X,t) = -Co2(X,t). 



(28) 



(29) 



(30) 



It is now clear that the quantities Cs^(X, r), which depend locally on time but nonlocally on space, play the role of 
effective fields because they encode the influence of the remaining particles on the particle located at X. They are 
also functionals of the scaled density 0oo(X, r), so that Eq. is self-consistent and can be solved numerically by 
iteration. Relation (|25p is a continuity equation for the particle density, whose current can be written as $i(X, r) — 

0ii(x,r)/r + 02i(x,T)/r2. 

Interestingly, for a = 1, Eq. f ([M]) ) can be recast to 



0ii(X,t) = -0oo(X,t)9 



'X 



^(X.r) 



+ V,.,t{X) 



where i^''°''[0oo] is the hard rod density functional of the instantaneous density 0oo(X, r). 

The term C02 vanishes in the limit a ^ and describes a tendency of the particles to form denser aggregates due 
to their inelasticity. 

Also notice that at the zero order in F^^, Eq. (P5)) may be rewritten as a dynamic density functional (DDF) 
equation^", since the only change with respect to the a = 1 caseJi is the presence of a prefactor (1 + a)/2 in the cqi 
term in Eq. (|26p . Hence, if we define an inelastic free energy density functional as Fa[p\ = Fideai[p\ + i^ + o:)/2Fi,xcess[p\ 
scaling the exact hard-rods excess of the equilibrium case, we would cast Eq. (|25p into a DDF equation for arbitrary 
values of a. A particular result would be that, always at the leading order in F^^, the equilibrium density profiles 
should be given by the minimum of such inelastic free energy density functional. The results in Figs. [5] and [7] are 
qualitatively consistent with that efi^ect since the reduction of the excess free energy reduces the oscillations. As we 
shall see below the equation of state for the uniform fluid (Eq. (|40)) ) predicts a lowering of the pressure with respect 
to the elastic case also consistent with such a scaling of the nonideal part of the free energy. 

Hereafter, we briefly derive some useful relations between the density proflle, the temperature and the pressure in 
the non uniform steady state where the current <i>i vanishes. 
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A. Steady state temperature profile 

In the limit t oo we determine the granular temperature profile, using the previous results and Eq. ([15 



T{X) = 1 



2 MX) 

rci>oa{x) 



2r0F 



g^iX, X + l)0oo(X + 1, r) + g^^X, X - l)</.oo(^ - 1) 



(31) 



where we suppressed the time argument and the functions of the single spatial argument have to be understood as 
their asymptotic limiting values when r oo. The constant 1 in the r.h.s. represents (in our reduced units) the 
heat-bath temperature, whereas the second term is the shift in the local temperature induced by collisions. In fact, 
it amounts to the product of three factors: the kinetic energy dissipated, the collision rate uje (Enskog collision 
frequency^'^) and the typical time of the heat-bath. The average Enskog frequency at each side of the particle 
located at X is: 



l^e{X±1) 

7 



<?2(X,X±1),/)00(X±1) 



(32) 



and in the case of a uniform system it reduces to the bulk Enskog frequency ue = ^vt p(y g2 1 ■ where vt is the 
thermal velocity of the gas. 



B. Steady state pressure profile. 



We now turn our attention to the pressure profile 11 (X, t), which can be separated into a kinetic and a coUisional 
contribution 



^{X)=^kUX) + ^cou{X) 



(33) 



The total pressure n(X, r) is implicitly determined from the momentum balance equation, obtained by considering 
Eq. (|14|) with = 1 in the t —* oo limit 



(34) 



(35) 



oX 

We identify the first term as 

nfc„(x) = y dVV^P{X, V) = $o(^) + 2$2(^) 

or using the results of the V expansion, we rewrite 

^Un{X) = (bm{X) + ^0i2(X) , = T{X)(boo{X) (36) 

where we used Eq. (|3ip to obtain the last equality. In Ref. fl3 we showed that the spatial derivative of the coUisional 
pressure is related to the collision integral via the relation 



dxTl,oii{X) = -Ci{X) = -coi(X) - ^cn(X) 



(37) 



By manipulating expressions ((26| and ((28)) (see Ref. [TO for details) we formally integrate Eq. ([37]) with the following 
result: 



^colliX) 



(l+«) r 



dzg-2,(X - (1 - z), X + z) X 



(38) 



|^Q„(X-(l-z))0oo(^ + ^) 

+ ^ \^MX - (1 - Z))(/)12(X + z)+ MX - (1 - z))MX + z) 

MX - (1 - z))MX + z)- MX - (1 - z))MX + z)\ } 
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In the case of constant density, the z-integration can be trivially performed, and we obtain 

(1 + a) (t^lo , 2(/)i2 



1 — 4>m V r Sm 



Finally, by using Eq. (|36p we cast the equation of state in the uniform non equilibrium steady state in the form 

(l + a) 000 



n = T0o 



1 



2 1-000 

Expressing the temperature as a function of the density (from Eq. (j3ip ) 

(1 

T =\ —^^52000 , 



(39) 



(40) 



(41) 



we see that Eq. (l40l) describes the lowering of the pressure due to the collisional reduction of the temperature and 
becomes the familiar hard-rod pressure equation for a — \. 



APPLICATIONS 



A. Temporal decay of a small density modulation 



We begin by considering the decay of an infinitesimal sinusoidal perturbation of wave-vector K with respect to a 
uniform density profile and how the relaxation time varies as a function of K. We assume that the sinusoidal density 
perturbation is small with respect to some uniform background density, $0 and write: 



(X,t) = $o + /5K(T)sin(A'X). 



(42) 



After some simple algebra we arrive at the following equation of evolution for the modulation: 



dpK{T) 
dr 



-R{K)pk{t) 



r 



-PK 



{r)[ 



2poCT- 



sin(i^) , 4(pofT)2 



K 



sin^{K/2) 



^ ^ 2il±^Po(T(l - cos(X)) 



[l-a^) (l + g)^ 
(l-«2) 



2 -^""2 



sin(ii') 4(no(T)^ . o,,^,„, 
2po(7 — ^ + sin^(i^/2) 



K 



r0F 



poo-(l + cos{K)) + (pqct) 



K 



1 



(1 + a) sui{K) 



-poa- 



K 



(43) 



with po = *I'o/(l — ^0) being the bulk hard rod reduced pressure. Taking the long wavelength limit iiT ^ 0, we obtain 
the following expression for the diffusion coefficient: 



D = hm = ± j 1 + (1 + a){p,a+pla^/2) 



l=^(2 + 3l±i^Po.)(po+P^./2)}. (44) 



As shown in Fig. [21 the relaxation time of the inelastic fluid, for small values of K , is longer than the corresponding 
quantity in the elastic fluid. This occurs because, being the granular temperature of the former lower, the diffusion 
is weaker. In addition, while the excluded volume favors diffusion with respect to the non interacting case, the 
inelasticity operates in the opposite direction. For K ^ Q the temperature field T(X, t) is anti-correlated with the 
density fluctuation p(X, r), and the local maxima of T correspond to the regions where the density is lower. On the 
other hand, for K larger and close to = tt, the temperature maxima occur in correspondence with the density 
maxima, thus we observe that the relaxation time of the inelastic system is shorter than the relaxation time for a = 1. 
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FIG. 2: Decay rate, R{K) of a small sinusoidal density perturbation around a constant background po = 0.68 as a function 
of the the reduced wave-vector K. The dimensionless friction constant is F = 5. The non monotonic dependence of R{K) on 
K increases as the density increases. Therefore collisions can accelerate or slow down the relaxation with respect to an ideal 
gas behavior. The dashed line represents the decay rate of a system of non interacting particles, the dash-dot line the Density 
functional result, the full line the system witha — 1, and the dot-line the inelastic system with a — 0.8. 



To validate the theory we shall compare the predictions of Eq. ([25|) with those obtained by a numerical solution of 
the dynamical equations for the trajectories of the particles, using the algorithm illustrated in Ref.llS- The comparison 
is performed by considering an ensemble of "noise" histories (typically 10'*) and averaging the observables over such an 
ensemble. The first category of checks concerns the homogeneous static properties of the system, namely temperature 
and pressure. The dependence of the temperature on the density predicted by Eq. (|4ip is shown in Fig. [Sj where it 
is displayed against the numerical results obtained at two different values of the coefficient of restitution. A similar 
comparison between pressure [Eq. (HO])] and the simulation results are reported in Fig. showing a satisfactory 
agreement. The second category concerns the inhomogeneous static properties, which can be probed by measuring 



FIG. 3: Uniform system: temperature versus reduced density for a — 0.9, a = 0.8 and F = 5. Comparison between Brownian 
dynamics simulations (points) and the predictions of our theory (lines). 

the response of the model to some specific non-uniform external perturbations. 

The third category of tests aims, instead, to probe some genuinely time-dependent properties of the system and 
we have chosen as examples the free expansion of a packet of particles initially localized in a narrow region and the 
escape of a packet from a potential well. 



B. Numerical tests 
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FIG. 4: Uniform system: pressure versus reduced density for a — 0.9, a = 0.9 and F = 5. Comparison between Brownian 
dynamics simulations (points) and the predictions of our theory (lines). 



C. Inhomogeneous Steady State Properties 



Soft repulsive potential. We consider, first, a fixed external potential of the form 

y(X) = Fotanh(X/C), (45) 

representing a soft repulsive wall located a,t X = and characterized by a softness parameter ^ = 0.2 and height 
Vb = 4. In Fig. [5] we compare the density profiles, obtained from the stationary solution of the dynamical Eq. (|25|). 
for two values of coefficient of restitution and for F = 5, with the corresponding profiles extracted from molecular 
dynamics simulation. The wall perturbs the fluid by inducing a non monotonic profile and a stationary state is 
achieved when hydrostatic equilibrium is reached. 

Both the MD and the integral equation reveal the same feature: near the wall the elastic density profile is slightly 
higher than the corresponding profile with a = 0.8. Intuitively such a difference can be understood by considering 
that the bulk pressures, approximately proportional to the corresponding values of the the density near the wall, 
display a similar difference. The temperature profile, shown in Fig. [6l varies non monotonically from the value inside 
the wall to a lower value in the bulk and is the signature of the non equilibrium nature of the system. 




FIG. 5: Density profiles in the presence of a repulsive soft wall ai x — 0, indicated by a dashed line, and implemented in the 
simulations through the potential V{X) — Vo tanh(X/^) with Vo = 4. Data refer to a = 1, a = 0.5 and F = 5. Points indicate 
the results of the simulations whilst the lines are the corresponding results from our theory. The agreement between simulation 
an theory is excellent and it is also interesting to note that the elastic system (a — 1) presents a higher peak near the wall 
corresponding to a larger pressure exerted, toward the wall, on each particle by the rest of the system. 

We study, now, the stationary profile induced by a static periodic external potential of the form 

V^(X) = Vocos(— X). (46) 
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FIG. 6: Temperature profiles in the presence of the same repulsive soft wall of Fig. [S]for a system with F = 5 and inelasticity 
a = 0.8, a = 0.5. 

For moderate values of the bulk packing fraction the asymptotic value of the induced density profile turns out to be 
modulated with the same period as the potential. The height of the peaks is lower than the corresponding height 
of the non-interacting case, because the hard-core repulsion tends to smear the particles over the wells away from 
the minimum energy configuration. On the other hand, one can appreciate a difference between the elastic and the 
inelastic case. The latter displays peaks slightly higher and narrower as shown in Fig. [71 a fingerprint of the tendency 
toward clustering induced by the inelasticity of collisions. 



0.55 




FIG. 7: Density profiles of a system with average density p = 0.45 in an external potential defined by Eq. (|46|l with w = 8 
and F = 5. Black circles correspond to simulation results with a = 0.8, while shaded circles refer to simulations with a = 1.0. 
Dashed and full lines indicate the corresponding results from the present theory for a = 1.0 and a — 0.8 respectively. 



D. Inhomogeneous Dynamical Properties 



We consider the free expansion of N hard-rods in the absence of external fields. In Fig. [51 we display the evolution 
of the variance of the positions of the particles, with respect to their center of mass, W = jjj^iii^i ~ Xcm)^) for 
different values of the inelasticity a. The average is meant over different and independent noise realizations. In the 
case of free particles D = F^^, we observe a linear growth of this quantity, which is well described by the diffusive law, 
W{t) = 2Dt. In agreement with our analytical prediction of Sec. VI, the coefficient D is a growing function of the 
coefficient of restitution and of the number of particles. The first effect can be interpreted by noticing that smaller 
values of a correspond to larger dissipation and thus to lower local kinetic temperatures. Since one expects D to be 
proportional to the ratio between temperature and friction coefficient a lowering of the kinetic temperature determines 
a decrease of the spreading. On the other hand, one can compare the spreading of the same initial configuration in 
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the case of non-interacting particles. Figure [8] shows that the ideal gas case corresponds to a diffusion slower of all 
cases where the hard core repulsion is at work. Finally, the importance of the corrections to the DDF equation, 
can be appreciated by rescaling the data corresponding to different values of F according to the formula W/T . The 
free particle case, of course gives a perfect collapse, whereas the interacting cases display increasing deviations as F 
decreases. 




FIG. 8: Free expansion of a cluster of 128 particles initially concentrated over a region of size 150cr. The curves represent the 
time growth of the variance of the particle distribution with respect to their center of mass for inelastic systems with different 
values of a but same F = 5. The inset shows the values of the diffusion coefficient estimated by the asymptotic slope of the 
curves. 




FIG. 9; Time decay of the number of particles initially located in a single well of width w for a = 0.8, F = 5. The same 
numerical experiment is repeated for different well sizes but maintaining fixed the barrier height. Broken lines represent the 
theoretical predictions and the symbols the simulations results obtained as an average over 300 independent runs. 

Finally, we consider the escape process of a group of particles from a well of the sinusoidal potential. This numerical 
experiment amounts to preparing initially a set of particles in a potential well and measuring the fraction nin(t) that 
populates the well at that instant. Figure [H reports the log-log plot of decay in time of ni„(t), obtained by averaging 
over 300 independent runs, for different well widths w. For comparison we also show the corresponding quantity 
computed through the DDFT (dashed lines) which correctly reproduces the dynamical features of the escape process. 
We see that the rate at which particles exit the original well decreases with w, because the particles spend more 
time in that well. The coUisionless particle systems generally displays a longer escape time than interacting systems, 
because the absence of excluded volume effects does not entail an effective reduction of the depth of the potential 
well, which, instead, is the relevant feature characterizing the escape experiments involving hard core particles. The 
inelasticity, on the other hand, tends to decrease the escape rate due to the energy loss caused by collisions, but this 
does not counterbalance the excluded volume effect so that the rate of the inelastic particles remains faster than the 
corresponding rate of the coUisionless model. 
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VI. CONCLUSIONS 

In this paper we have derived a method to study the dynamics of an assembly of particles interacting inelastically, 
and driven by a stochastic thermostat. We have found that the particles adopt spatial configurations which are very 
close to those of an equilibrium system, in spite of the fact that our system is driven and dissipative. The reason for 
such a similarity is twofold as suggested by a recent study of Reis et al4^: the homogeneous energy feeding mechanism 
and the importance of the repulsive forces. However, the present work shows that there is no need to invoke entropic 
forces to explain the observed inhomogeneities. A kinetic approach, in which the short range repulsion is suitably 
accounted for by means of a suitable treatment of spatial correlations, predicts fairly well the observed structural 
properties^. Our theory indicates that the steady state configurations occur not as a result of the minimization of 
some hypothetical coarse grained free energy functional but as a result of the competition between the uniform energy 
injection and the energy dissipation. These two effects are described by an Enskog collision operator and by a a 
Fokker-Planck collision operator, respectively. The evolution of the phase space distribution function is thus governed 
by a Fokker-Planck-Enskog (FPE) equation, whose solution still remains an extremely difficult task for dense fluids 
due to the complexity of the collision kernel and to the computer resources needed to resolve the distribution function. 
However, when the friction is sufficiently high one can derive a simpler description by an iterative elimination of the 
fast degrees of freedom, such as the velocities of the particles. Such a procedure is based on the intuition that these 
achieve locally their equilibrium distribution, whereas the positions evolve more slowly. Truncating this iteration at 
the first order in the inverse friction parameter is equivalent to approximating the velocity distribution functions 
by Maxwellians at temperatures equal to that of the heat bath. Further terms, associated with non Maxwellian 
contributions to the velocity distribution function, are included in the expansion and contribute to the evolution. The 
result is a self-consistent time dependent equation for the local density, where the "internal field" is determined by 
the density itself either through the standard hard-rod entropic contribution or by the velocity and energy currents 
generated by spatial density gradients. 

At a technical level Eq. ([^5]) is derived by applying a multiple time-scale method to the Fokker-Planck-Enskog 
equation. The resulting equation bears a strong similarity with the Dynamic Density Functional equation, but it is 
not based on the notion of coarse grained Free Energy, a concept which cannot be applied to open non equilibrium 
systems, such as the inelastic fluid we have studied in this work. Equation ([^5)1 is a density functional equation for 
the density field (f>oQ{X,T) where its evolution depends on a functional of (/>oo itself and its derivatives. It is local in 
time, but does not possess a generating Liapunov functional, therefore, we cannot prove that the associated dynamics 
minimizes some cost function. 

In more detail we have found that: 

a) the inelasticity induces changes even in the stationary properties of the fluid with respect to the elastic reference 
system to order F"^. The changes can be observed both in the velocity distribution and in the structure of the non 
uniform fluid. 

b) Whereas in the description of a colloidal fluid the friction F is relevant only for the relaxation properties, in the 
inelastic fluid F determines its stationary properties also. 

c) The theory holds in the region useIi < when the typical time-scale of the heat bath is shorter than the Enskog 
collision frequency. 

As far as future perspectives are concerned the method can be generalized to higher dimensions and different types of 
inter-particle forces and to systems with a non uniform distribution of heat sources^. A second type of generalization 
consists in performing the same multiscale expansion at the level of the two particle phase-space distribution function, 
by truncating the BBGKY hierarchy one step further, which would allow us to compute self-consistently the pair 
correlation function of the system. 
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APPENDIX A: COLLISION INTEGRALS 

In this appendix we show how to perform the velocity integrations and reduce the collision integrals to simple 
functions of space and time only. Using the definition of collision integral given in the text by Eq. ([5]) and setting 
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u = V2 — V, we obtain the following explicit expression: 

CniX,T) ^ g2{X,X + l)[f dV^ln{V)[^ duuP {X , V, t)P {X + 1 , u + V, t) 

+ X; ( duuP{X,V + su,t)P{X + l,V + qu,T)\\ 
Jo J J 

1 r 



32(^,^-1) / dV 



duuP{X, V + SM, T)PiX -l,V + qu,T) 



duuP{X,V,T)P{X -1,u + V,t) I 



(Al) 



where s = {l + a)/2a, q = —(1 — a)/2a, fio{V) — 1, fJ.i{V) = V and fJ.2{V) — After substituting the expression 

of P{X, V, t) in terms of its partial amplitudes into Eq. (jAl|) . one can eliminate the velocities obtaining: 



C„(X,r) = g2{X,X + l)Y, 
- g2{X,X-l)Y, 



Af(^)(a = l) + i.7V(")(a) 



$^(X,T)$„(X + l,r) 
<i>^(X,T)$,(X-l,r). 



where the matrix elements Mu"-* and are defined as: 



duu / dT/yLt„(F)iJ^(F + gu)iJi.(F + su) 



(A2) 



(A3) 



M^rJia)^ duu dVfin{V)Hf,{V + qu)H,{V + su). 

Jo J -oo 

and have the symmetry property 

so that it is sufficient to calculate only the matrix elements of M^^i) {a) in order to compute Eq. (jA2 



A4l'(a) 



-a/2 
2^(1 + 2") 



l~2a 



a/2 



\{l-a) 




-^(l-2a2) -i(l + 2a-3a2) ^(l-Sa + Ga^) 
i(l-2a-3Q2) _3^(i_2a2) i(l-a)(l + 3a) 
_^(l + 8a + 6a2) -1(1 + «)(!- 3«) -^(1 - 2^2) 



In addition, we verify that in the case rt = the following combinations vanish: 



(A4) 



(A5) 



so that Co{X,t) — 0, since collisions conserve the number of particles. 
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APPENDIX B: MULTIPLE TIME SCALE METHOD. 



Our previous work'^-'- has extended to the case of cohiding particles a method to derive the Smoluchowski equation 
starting from the Kramers equation. It was originaUy proposed in the 1970s for a gas of non-interacting particles by 
Titulaer"'^ and nicely reviewed by Bocquet^i^-. 

It represents a particular application of multiple time-scale analysis^ designed to handle singular perturbations. 
In the present case the singularity stems from the fact that when F > > 1 the time derivative occurs among the small 
terms of Eq. (fT4|) . 

Because the inelasticity brings about some remarkable new features we shall report the derivation of the salient 
parts of the multiple scale method in this particular case. The multiple time-scale analysis introduces a set of auxiliary 
time scales t„ — r""T, with n = 0, 1,2, .... The t„ are treated as independent variables so that the time derivative 
with respect to r is replaced by 



d _ d 19 15 
97 ^ 9^^f 9^^f2 9^ 



(Bl) 



The partial amplitudes 4>^(X, r) and the collision terms C^(X, r) are also treated as functions of the auxiliary time 
scales and expanded perturbatively as: 



$i.(X,ro,ri,T2, ..) = ^ — 0„^(X,To,Tl,T2, ..) 



Tl = 



and 



OO 

'^(X, To, Ti, T2, ..) = ^ — C„y(X, ro,Ti,T2, 



(B2) 



(B3) 



n=0 



By substituting Eqs. (jBip - (|B3p into Eq. and equating equal powers of F, one obtains iteratively a series of 
equations which must be satisfied by the coefhcients <^sv(X^ tq,ti,T2, ..) and Csu{X, tq,ti,T2, ..). The latter coefficients 
are obtained using the formula: 



l+m=s lJ.,v 



N\rJ{a = l) + ^MirXa) 



52(X,X-1) 



(B4) 



Notice that the c^i/'s are functionals of the ^si/'s. 
We begin with the order F": 



(B5) 



having the solution (j)^^ = for v ^ which inserted in Eq. ^ determines the expansion coefficients of order F ^ in 
terms of (f>Q^ = 0. 



L 



FP 



5<Aoo 



'I2H2 + 0l3i?3 + 014-^4 + 

^^0 + Dx(l>ooHi — coi-ffi — C02H2 — C03H3 — C04H4. 



(B6) 



where we have employed the abbreviation Dx = {dx — F{X)). We also perform our expansion by setting (pso = for 
all s > 0. By equating the coefficients of the same H^, in Eq. (jB7P we find the following relations 







(B7) 



hi = -Dx4>m + coi 



(B8) 
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and for v > 1 

<j)iu = -^co^. (B9) 
The procedure can be iterated to the order F"^, writing the equation 

LFp[Y,(t^2uH,] -^^i/. + ^i/o + ^A-01li^2 + 9x011^0+ (BIO) 

L»js:0i2i?3 + 29x012^^1 + Dx(l)i3H4 + 39x013^^2 + Dx(t>iiH^ + 49x0i4i?3 - ^ ci^H^ 

v>\ 

which leads to the foUowing conditions 



-9x011 (Bll) 



9to 

9012 



" -021 - 29x012 + cii =0 (B12) 



9to 

and 

9013 

9to 



= -2022 --0x011- 39x013 + ci2 = 0, (B13) 



-3023 --0x012-49x014+013 = 0. (B14) 



The amphtude 0ii, being a functional of 0oo, does not depend on tq. Hence, the l.h.s. of Eq. (|B12p vanishes and we 
find 



Similarly, we obtain for v > \: 



Explicitly we write 



21 =-29x012 +C11. (B15) 



hu = --[-0x01(1.-1) + (j^ + l)9x0i(,.+i) -I- cxv\ (B16) 



^22 = -^[-Ox -0x000 + -Dxcol -I- 9xCo3 - C12] (B17) 

023 = -^[^-OxC02 + 9xC04 - C13] (B18) 

024 = -^[^-OxC03 + 9xC05 - Cm]. (B19) 



Finally, by equating the coefhcients of _ffo, we arrive to the equation: 

9000 



9ti 



= 9x[-Dx0oo-coi]. (B20) 



In order to carry out the derivative of 0oo with respect to the time r2, we iterate the procedure to the order F ^ 
by writing: 

^J^p f V 03.i/J = 7^ E + IT- E -^i-^- + (B21) 

L^i ^ ^-"^t^i ^^1^1 

+ E['^^'^2i/^fi.+l + i^9x02<.-H^,y-l] - E '^2!/ 
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and equating the coefRcients of Ho{V) we obtain the following equation 

^ = ~dx<p2i = -dxcii + aico2 (B22) 

where the second equality follows from eqs. l|B9[) and (jBlsp . 

We, now, collect together the different orders in r„, given by Eqs. (|Bip . (jBlip and (|B22[) . thus restoring the original 
physical time r to find the evolution equation for the density amplitude (1251) . 
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